# Indicators
allIndicators <- read_csv(here::here("Processed_Indicators/allIndicators.csv"))
allIndicators$Season <- factor(allIndicators$Season, levels= c("spring", "summer", "fall", "winter", "all"))
allIndicators$stat_area <- factor(allIndicators$stat_area, levels= c("511", "512", "513", "511-513"))
indicators <- unique(allIndicators$stat_area)
plot_fun <- function(x, indicator){
  x %>% 
    filter(Indicator == indicator) %>% 
    ggplot() + 
    geom_line(aes(Year, Value, col = stat_area)) +
    theme_bw() +
    facet_wrap(~Season)
}

sumLms_fun <- function(indicator, season, stat_areas){
  df <- allIndicators %>% 
    filter(Indicator == indicator,
           Season == season,
           stat_area == stat_areas)
  
  lm1 <- lm(Value ~ Year, data = df)
  return(lm1)
}


acf_plot <- function(indicator, season, stat_areas){
    df <- allIndicators %>% 
    filter(Indicator == indicator,
           Season == season,
           stat_area == stat_areas) %>% 
      select(Year, Value)
  acf(df)
}

pscore_fun <- function(indicator, season, stat_areas){
  df <- allIndicators %>% 
    filter(Indicator == indicator,
           Season == season,
           stat_area == stat_areas)
  
  lm1 <- lm(Value ~ Year, data = df)
  pscore <- segmented::pscore.test(lm1)
  return(pscore)
}

davies_fun <- function(indicator, season, stat_areas){
  df <- allIndicators %>% 
    filter(Indicator == indicator,
           Season == season,
           stat_area == stat_areas)
  
  lm1 <- lm(Value ~ Year, data = df)
  davies <- segmented::davies.test(lm1)
  return(davies)
}

slope_bp_fun <- function(indicator, season, stat_areas, Npsi){
    df <- allIndicators %>% 
    filter(Indicator == indicator,
           Season == season,
           stat_area == stat_areas)
  
  lm1 <- lm(Value ~ Year, data = df)
  
  segmented::segmented(lm1, seg.Z = ~Year, npsi = Npsi)
}

mean_bp_fun <- function(indicator, season, stat_areas, models){
    df <- allIndicators %>% 
    filter(Indicator == indicator,
           Season == season,
           stat_area == stat_areas)
  
  lm1 <- lm(Value ~ Year, data = df)
  temp_mcp1 <- mcp::mcp(models, data = df, par_x = "Year")
  return(temp_mcp1)
}

Developing Indicators of Habitat and Ecosystem Change in the Gulf of Maine

This report includes the initial analysis of the ecosystem and habitat indicators. Linear regression was used to assess trends over the length of the time series. A trend breakpoint analysis and a mean breakpoint analysis were run on the indicators. A breakpoint found in the slope of the line indicates a change in the trend of the data. A breakpoint in the mean of the date may indicate a regime shift or change of the overall state of the system. These analyses were run seasonally within statistical areas, when available, and as yearly averages across the entire study domain.

Indicators

The following indicators are used in this report

  • Surface (OISSTv2.1, FVCOM NECOFS)
  • Bottom temperature anomalies (FVCOM NECOFS)

OISST Temperature

Time series plots

plot_fun(allIndicators, "oisst")

Linear models

lms <- sumLms_fun("oisst", "all", "511-513")
plot(lms)

summary(lms)
## 
## Call:
## lm(formula = Value ~ Year, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.15193 -0.32238  0.06949  0.35912  1.38051 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -78.585742   9.579605  -8.203 3.80e-12 ***
## Year          0.039408   0.004789   8.230 3.38e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4944 on 78 degrees of freedom
## Multiple R-squared:  0.4648, Adjusted R-squared:  0.4579 
## F-statistic: 67.73 on 1 and 78 DF,  p-value: 3.376e-12
acf_plot("oisst", "all", "511-513")

Slope breakpoint tests

The pscore test and davies tests indicate whether a breakpoint in the slope exists. The davies test is more sensitive to data that follow a more sinusoidal curve but less sensitive to linear changes. A significant p-value indicates there is a breakpoint in the data. A non-significant p-value indicates a breakpoint is not present.

pscore_fun("oisst", "all", "511-513")
## 
##  Score test for one/two changes in the slope
## 
## data:  formula = Value ~ Year 
## breakpoint for variable = Year 
## model = gaussian , link = identity , method = lm
## observed value = 3.3271, n.points = 10, p-value = 0.00134
## alternative hypothesis: two.sided   (1 breakpoint)
davies_fun("oisst", "all", "511-513")
## 
##  Davies' test for a change in the slope
## 
## data:  formula = Value ~ Year ,   method = lm 
## model = gaussian , link = identity  
## segmented variable = Year
## 'best' at = 1994, n.points = 8, p-value = 0.0005195
## alternative hypothesis: two.sided

Slope breakpoint

bp <- slope_bp_fun("oisst", "all", "511-513", 1)
summary(bp)
## 
##  ***Regression Model with Segmented Relationship(s)***
## 
## Call: 
## segmented.lm(obj = lm1, seg.Z = ~Year, npsi = Npsi)
## 
## Estimated Break-Point(s):
##            Est. St.Err
## psi1.Year 1992  2.029
## 
## Meaningful coefficients of the linear terms:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 86.40086   52.38220   1.649    0.103
## Year        -0.04359    0.02637  -1.653    0.102
## U1.Year      0.10377    0.02738   3.790       NA
## 
## Residual standard error: 0.4459 on 76 degrees of freedom
## Multiple R-Squared: 0.5757,  Adjusted R-squared: 0.559 
## 
## Convergence attained in 1 iter. (rel. change 3.4373e-08)
plot(bp)
df <- allIndicators %>% 
  filter(Indicator == "oisst",
         Season == "all",
         stat_area == "511-513")
points(x = df$Year, y = df$Value)

Mean breakpoint

mean_cp <- mean_bp_fun("oisst", "all", "511-513", list(Value~1, 1~1))
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 80
##    Unobserved stochastic nodes: 4
##    Total graph size: 580
## 
## Initializing model
summary(mean_cp)
## Family: gaussian(link = 'identity')
## Iterations: 9000 from 3 chains.
## Segments:
##   1: Value ~ 1
##   2: Value ~ 1 ~ 1
## 
## Population-level parameters:
##     name     mean   lower   upper Rhat n.eff
##     cp_1 2009.354 2008.53 2.0e+03  1.2  1198
##    int_1   -0.073   -0.19 4.1e-02  1.0  4238
##    int_2    1.089    0.91 1.3e+00  1.0  2597
##  sigma_1    0.427    0.36 5.0e-01  1.0  4734
plot(mean_cp)

FVCOM SST Temperature

Time series plots

plot_fun(allIndicators, "fvcom_sst")

Linear models

lms <- sumLms_fun("fvcom_sst", "all", "511-513")
plot(lms)

summary(lms)
## 
## Call:
## lm(formula = Value ~ Year, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.35607 -0.35722 -0.03152  0.43646  1.60338 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -64.956743  10.689314  -6.077 3.94e-08 ***
## Year          0.032429   0.005345   6.068 4.10e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5726 on 80 degrees of freedom
## Multiple R-squared:  0.3152, Adjusted R-squared:  0.3066 
## F-statistic: 36.82 on 1 and 80 DF,  p-value: 4.097e-08
acf_plot("fvcom_sst", "all", "511-513")

Slope breakpoint Tests

The pscore test and davies tests indicate whether a breakpoint in the slope exists. The davies test is more sensitive to data that follow a more sinusoidal curve but less sensitive to linear changes. A significant p-value indicates there is a breakpoint in the data. A non-significant p-value indicates a breakpoint is not present.

pscore_fun("fvcom_sst", "all", "511-513")
## 
##  Score test for one/two changes in the slope
## 
## data:  formula = Value ~ Year 
## breakpoint for variable = Year 
## model = gaussian , link = identity , method = lm
## observed value = 3.7145, n.points = 10, p-value = 0.0003754
## alternative hypothesis: two.sided   (1 breakpoint)
davies_fun("fvcom_sst", "all", "511-513")
## 
##  Davies' test for a change in the slope
## 
## data:  formula = Value ~ Year ,   method = lm 
## model = gaussian , link = identity  
## segmented variable = Year
## 'best' at = 1988.9, n.points = 8, p-value = 0.0003583
## alternative hypothesis: two.sided

Slope breakpoint

bp <- slope_bp_fun("fvcom_sst", "all", "511-513", 1)
summary(bp)
## 
##  ***Regression Model with Segmented Relationship(s)***
## 
## Call: 
## segmented.lm(obj = lm1, seg.Z = ~Year, npsi = Npsi)
## 
## Estimated Break-Point(s):
##                Est. St.Err
## psi1.Year 1989.496  1.977
## 
## Meaningful coefficients of the linear terms:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 163.65082   80.13934   2.042   0.0445 *
## Year         -0.08267    0.04038  -2.047   0.0440 *
## U1.Year       0.13529    0.04105   3.296       NA  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5187 on 78 degrees of freedom
## Multiple R-Squared: 0.4521,  Adjusted R-squared: 0.431 
## 
## Convergence attained in 2 iter. (rel. change 0)
plot(bp)
df <- allIndicators %>% 
  filter(Indicator == "fvcom_sst",
         Season == "all",
         stat_area == "511-513")
points(x = df$Year, y = df$Value)

Mean breakpoint

mean_cp <- mean_bp_fun("fvcom_sst", "all", "511-513", list(Value~1, 1~1))
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 82
##    Unobserved stochastic nodes: 4
##    Total graph size: 594
## 
## Initializing model
summary(mean_cp)
## Family: gaussian(link = 'identity')
## Iterations: 9000 from 3 chains.
## Segments:
##   1: Value ~ 1
##   2: Value ~ 1 ~ 1
## 
## Population-level parameters:
##     name    mean   lower   upper Rhat n.eff
##     cp_1 2009.47 2009.02 2010.00    1  4233
##    int_1   -0.43   -0.54   -0.31    1  5679
##    int_2    0.79    0.61    0.97    1  5404
##  sigma_1    0.43    0.37    0.50    1  5211
plot(mean_cp)

FVCOM BT Temperature

Time series plots

plot_fun(allIndicators, "fvcom_bt")

Linear models

lms <- sumLms_fun("fvcom_bt", "all", "511-513")
plot(lms)

summary(lms)
## 
## Call:
## lm(formula = Value ~ Year, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.95338 -0.35913  0.03703  0.33016  1.27851 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -27.638215   8.590416  -3.217  0.00187 **
## Year          0.013804   0.004295   3.214  0.00189 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4602 on 80 degrees of freedom
## Multiple R-squared:  0.1144, Adjusted R-squared:  0.1033 
## F-statistic: 10.33 on 1 and 80 DF,  p-value: 0.001888
acf_plot("fvcom_bt", "all", "511-513")

Slope breakpoint Tests

The pscore test and davies tests indicate whether a breakpoint in the slope exists. The davies test is more sensitive to data that follow a more sinusoidal curve but less sensitive to linear changes. A significant p-value indicates there is a breakpoint in the data. A non-significant p-value indicates a breakpoint is not present.

pscore_fun("fvcom_bt", "all", "511-513")
## 
##  Score test for one/two changes in the slope
## 
## data:  formula = Value ~ Year 
## breakpoint for variable = Year 
## model = gaussian , link = identity , method = lm
## observed value = 2.7247, n.points = 10, p-value = 0.007903
## alternative hypothesis: two.sided   (1 breakpoint)
davies_fun("fvcom_bt", "all", "511-513")
## 
##  Davies' test for a change in the slope
## 
## data:  formula = Value ~ Year ,   method = lm 
## model = gaussian , link = identity  
## segmented variable = Year
## 'best' at = 1993.3, n.points = 8, p-value = 0.0244
## alternative hypothesis: two.sided

Slope breakpoint

bp <- slope_bp_fun("fvcom_bt", "all", "511-513", 1)
summary(bp)
## 
##  ***Regression Model with Segmented Relationship(s)***
## 
## Call: 
## segmented.lm(obj = lm1, seg.Z = ~Year, npsi = Npsi)
## 
## Estimated Break-Point(s):
##            Est. St.Err
## psi1.Year 2004  3.813
## 
## Meaningful coefficients of the linear terms:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.064436  17.203649   0.585    0.560
## Year        -0.005135   0.008636  -0.595    0.554
## U1.Year      0.053291   0.018968   2.810       NA
## 
## Residual standard error: 0.4404 on 78 degrees of freedom
## Multiple R-Squared: 0.2093,  Adjusted R-squared: 0.1789 
## 
## Convergence attained in 1 iter. (rel. change -1.3997e-08)
plot(bp)
df <- allIndicators %>% 
  filter(Indicator == "fvcom_bt",
         Season == "all",
         stat_area == "511-513")
points(x = df$Year, y = df$Value)

Mean breakpoint

mean_cp <- mean_bp_fun("fvcom_bt", "all", "511-513", list(Value~1, 1~1))
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 82
##    Unobserved stochastic nodes: 4
##    Total graph size: 594
## 
## Initializing model
summary(mean_cp)
## Family: gaussian(link = 'identity')
## Iterations: 9000 from 3 chains.
## Segments:
##   1: Value ~ 1
##   2: Value ~ 1 ~ 1
## 
## Population-level parameters:
##     name    mean   lower    upper Rhat n.eff
##     cp_1 2009.22 2008.00 2010.187    1  3279
##    int_1   -0.19   -0.30   -0.085    1  5191
##    int_2    0.40    0.22    0.573    1  4080
##  sigma_1    0.41    0.35    0.479    1  5006
plot(mean_cp)